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CROSS REFERENCE TO RELATED APPLICATIONS 



[001] This application claims priority from United States Provisional App. Ser. No. 
60/416,140 filed on October 4, 2002. 

5 

FIELD OF THE INVENTION 



[002] This invention relates to the field of seismic data processing and, more 
particularly, to the use migration for determining subsurface earth structure represented 
10 by a 3-D volume of data for identifying structural and stratigraphic features in three 
dimensions. 



BACKGROUND OF THE INVENTION 



15 [003] Searching for subsurface mineral and hydrocarbon deposits comprises data 

acquisition, analysis, and interpretation procedures. Data acquisition involves energy 
sources generating signals propagating into the earth and reflecting "from subsurface 
geologic structures. The signals received are recorded by receivers on or near the surface 
of the earth. The received signals are stored as time series (seismic traces) that consist of 

20 amplitudes of acoustic energy which vary as a function of time, receiver position, and 

source position and, most importantly, vary as a function of the physical properties of the 
structures from which the signals reflect. The data are generally processed to create 
volumes of acoustic images from which data analysts (interpreters) create maps and 
images of the subsurface of the earth. 
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[004] Data processing involves procedures that vary depending on the nature of the data 
acquired and the geological structure being investigated. A typical seismic data 
processing effort produces images of geologic structure. The final product of data 
5 processing sequence depends on the accuracy of these analysis procedures. 

[005] Processed seismic data are interpreted to make maps of subsurface geologic 
structure to aid decisions for subsurface mineral exploration. The interpreter's task is to 
assess the likelihood that subsurface hydrocarbon deposits are present. The assessment 
10 will lead to an understanding of the regional subsurface geology, important main 

structural features, faults, synclines and anticlines. Maps and models of the subsurface, 
both in 2D and 3D representations are developed from the seismic data interpretations. 
As is well known in the art, the quality and accuracy of the seismic data processing has a 
significant impact on the accuracy and usefulness of the interpreted data. 

15 

[006] High quality data processing greatly simplifies data interpretation, since resources 
can be focused on the geologic structure since subsurface imaging can be made less 
ambiguous. Unfortunately, three dimensional geophysical data processing and/or 
modeling frequently require large computation expenses, and practitioners are forced to 
20 simplify the data processing effort as much as possible to reduce analysis time and cost. 

[007] The sheer volume of data impacts data processing considerations. Seismic survey 
data sets can involve hundreds of thousands of source locations, with each source 
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location associated with many hundreds more receiver locations. Each input/output data 
transfer demand burdens resources independent of the computation burden. 

[008] There have been several different approaches to manage these computational 
5 resource burdens. These approaches relate to the manner in which the data acquisition 
exercise is designed and carried out, as well as to assumptions made during data 
processing. The us 3 of available a priori geologic and geophysical information can 
facilitate the minimization of the seismic data acquisition effort. Such a minimization of 
resources reduces the amount of data that is acquired by reducing the acquisition effort. 

10 

[009] Minimization of the computational effort is often implemented during data 
processing. Compromises often required during data acquisition and processing can result 
in ambiguous and/or inaccurate subsurface images. Because little is generally known of 
the geologic structure being investigated, the interpreter will not know the extent the 
1 5 images are erroneous. 

[0010] It is not uncommon for significant computer resources to be involved when large 
or complex data volumes are processed, often involving weeks or months of actual 
computer processing time. The recent availability of massively parallel processor 
20 computers offers a significant opportunity to reduce overall processing times. Massively 
parallel processors (MPPs) can have multiple central processing units (CPUs) which can 
perform simultaneous computations. By efficient use of these CPUs, projects that took 
weeks or months of resource time previously can be reduced to a few days or a few 
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hours. These advantages can be enhanced further when efficient algorithms are included 
in the MPP software. 

[0011] Computational algorithms have previously been written for prior seismic analysis 
5 routines using single or just a few processors, usually using sequential computing. 

Sequential computing performs single procedures at any given time. Options for 
obtaining enhanced performance are limited when few processors are available. 

[0012] MPP computing machines offer an obvious computation advantages. The total 
10 time required to process a dataset can be reduced by dividing the work to be done among 
the various CPUs or CPU clusters in manner such that each CPU performs useful work 
while other CPUs also work in parallel. 

[0013] Seismic data consists generally of a large number of individual traces, each 
15 recorded somewhat independently of the other traces. Logically enough, sequential 

computing methods that require the processing focus to be placed on a single calculation 
at a time adapt well to analysis of these independent traces using parallel processing. This 
is true even though computational bottlenecks may exist. For example, portions of the 
processing sequence may require relatively more computation time than other portions, 
20 must be completed before other calculations may proceed, or may rely on input data other 
than seismic data, *br example traveltimes. Even though parallel processing significantly 
decreases total processing times, efficient processing algorithms can also contribute 
significant time reductions. 
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[0014] For large datasets frequency domain methods offer significant computation cost 
savings. An option for overcoming increased computational loads is to employ efficient 
frequency domain migration algorithms. 

[0015] It would therefore be desirable to have a system and method that is able to 
provide a migration in the frequency-wavenumber domain on parallel computer in a cost- 
effective manner. The present invention satisfies this need. 

SUMMARY OF THE INVENTION 

[0016] The method and system of the present inventions provides a new approach for 
selecting frequencies comprising a minimal or limited frequency data set for imaging 
leading to an order of magnitude decrease in computation resources. Seismic data are 
acquired and transformed a frequency domain. The data may be gathered in any manner. 
A plurality of minimal frequency subsets are formed comprising selected frequencies. 
The frequencies may be selected so that more frequencies are selected for shallower 
image positions. The selected frequencies are imaged to form frequency images. The 
selected frequencies may be individually weighted assure a full and balanced bandwidth. 
The frequency images are composited to form an intermediate or final composited 
migration image. 
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BRIEF DESCRIPTION OF THE DRAWINGS 



[0017] The present invention and its advantages will be better understood by referring to 
the following detailed description and the attached drawings in which: 

Figure 1: Illustrates wavenumbers imaged for a given frequency. 

Figure 2: Illustrates a full-bandwidth shot-domain migration. 

Figure 3: Illustrates two-dimensional shot-domain migration. 

Figure 4: Illustrates two-dimensional shot-domain migration. 

Figure 5: Illustrates two-dimensional shot-domain migration. 

Figure 6: illustrates variable downward continuation with intermediate spline 

interpolation. 

Figure 7: illustrates a depth slice from a common-azimuth-limited-frequency migration 
of a large 3D Gulf of Mexico data set. 

Figure 8: illustrates a flow chart of an embodiment of the present invention. 

Figure 9: illustrates a computer system for carrying out an embodiment of the invention. 

[0018] While the invention will be described in connection with its preferred 
embodiments, it will be understood that the invention is not limited thereto. On the 
contrary, it is intended to cover all alternatives, modifications, and equivalents which 
may be included within the spirit and scope of the invention, as defined by the appended 
claims. 
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DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 



[0019] Zero-offset-migration/inversion of single frequency data produces a subsurface 
image with spatial wavenumbers confined to an Ewald sphere defined by a simple 
dispersion relation. For non-zero-offset, the Ewald sphere becomes an Ewald doughnut. 
A single frequency images a wavenumber volume rather than just a spherical shell. This 
fact can be exploited to substantially reduce the number of frequencies needed to 
accurately map subsurface reflectors. This invention provides a method and apparatus for 
calculating a sparse set of frequencies for use in frequency-slice migration algorithms. 
The migration process remains the same, but the smaller set of frequencies results in an 
algorithm an order of magnitude faster than more traditional implementations. A small 
collection of case studies compares the limited frequency approach to traditional 
methods. Finally, we discuss various strategies for parallel implementation in distributed 
computational environments. 

[0020] The standard approach to frequency slice migration algorithms (phaseshift, 
split-step, ({fyx) finite difference, phase-screen) is to separately migrate each frequency 
independently, using an increment df defined by at least twice the total input time series 
length. If the trace length is 7, then wrap-around is avoided by simply computing the FFT 

over the length IT. This results in a frequency increment of df = ^ and means that the 

number of frequencies that must be migrated is 
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H fre q = 



df 



(1) 



For the average data set, n freq is typically one-half to one-third of the number of time 
samples in each trace, and ranges from a few hundred to around one-thousand. 

5 

[0021] Recently, Sirgue, Pratt, Plessix and Mulder (Sirgue and Pratt (a), 2001; Sirgue and 
Pratt (b), 2001; Plessix et. al., 2001) propose that full waveform inversion can be 
accomplished using an extremely limited discrete set of frequencies. In one-dimension, 
they show that the alternative set of frequencies can be defined by letting be the 

10 maximum depth of the migration, the maximum offset in the data, and 



1 



Then, in a homogeneous medium with velocity, v, the migration of each frequency results 
15 in an image with a spatial bandwidth of 

V V 



20 



Note that slower velocity media produce wider spatial bandwidths per frequency thai . 
faster ones, but the lower end of the spatial bandwidth is controlled by a. 
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[0022] It is natural to ask if this idea might be effective in the simpler case of downward 
continuation based wave-equation migration. More importantly, it is of interest to 
determine the extent to which this concept applies in higher dimensions. Extension of the 
methodology to higher dimensions can be accomplished in several ways. An immediate 
and positive conclusion is based on the so called Ewald sphere of Clayton and Stolt 
(1981). For example, in a 3D zero-offset setting, migration or inversion of a single 
frequency resolves the spatial wavenumbers lying in a spherical dispersion- relation- 
defined shell. In the nonzero offset case the shell becomes a volume, and the associated 
dispersion-relation-controlled figure might be described as an Ewald doughnut. The 
definition of a is slightly different but the fact that a volume of wavenumbers is imaged 
by a single frequency is clearly valid. 



[0023] A more detailed explanation can be based on the wave- number-frequency- 
domain-double-square-root equation. After mathematical rationalization, one has 



k 2 - 



-) = ft 2 2 + ft If (1 + tan 2 /) (4) 



where k 0 = — is the temporal wavenumber, k 2 the vertical wavenumber, k. the lateral 
v 

k 

spatial wavenumber vector, k h the offset wavenumber vector, and y = tan 1 (— ) is the 

K 

opening angle defined by a pair of rays impinging and reflecting from a subsurface 
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reflector. Note that for any fixed angular frequency <y , this quartic defines a five 
dimensional volume. The resulting figure is definitely not a shell. What is even more 
important is the fact that for the usual seismic acquisition geometry, opening angles y 
are large for shallow events but small for deeper events. As a rule of thumb, y , can be 
5 assumed to be zero when the depth of penetration is equal to or greater than the 

maximum offset in the acquisition scheme. 

[0024] To get an idea of the shape of the wave-number solid, it is convenient to consider 
the projection-slice theorem, 

10 

p(k s9 k r9 co)=k 2 Q o(k h \ (5) 

where P is the downward continued wavefield, k s and k r are the source and receiver 

wavenumbers, and O is the final output image. In two dimensions, Figure 1 provides a 
15 rough idea of the range of wavenumbers that are imaged by a given fixed frequency. The 

area of the large semicircle 101 indicates imaged wavenumbers. The smaller semicircles 
103 represent unimaged wavenumbers. As the frequency increases, the areal extent of the 
wavenumber images increase geometrically. The bad news is that higher frequencies are 
not weighted in the correct manner. This issue must be addressed in order to maintain the 
20 overall quality of the final image. The actual number of frequencies required to produce a 
high quality image decreases as a function of penetration depth. Another way to say this 
is that any limited frequency strategy should use more frequencies shallow than are used 
deep. 

COR- 1058 11 



[0025] For a given bandwidth (/ 0 »/i) choose 



N = 



log f 0 -log/, 
loga(z) 



(6) 



where 



a(z) = 



1 1 



1 + 



fO 2 Vl + tan 2 r (z) 



is chosen to vary with z. The number of frequencies to be migrated will then vary 
vertically and begin with the full bandwidth and end with a small percentage of the 
original set. 

[0026] Start by migrating only the frequencies 



from n = 0 to n = N . As the downward continuation proceeds, decrease N, through 
a(z) so as to reduce the number of frequencies required for the migration. As the minimal 
frequency subset changes it may be necessary to use interpolation to continue the process 
in an accurate manner, but the result is a very computationally efficient algorithm. At 
each depth level, individual frequencies can be weighted to assure that a full and 
balanced bandwidth is imaged. The result is a set of minimal frequency data subsets that 



/„=«(*)■"/. 



(7) 



COR-1058 



12 



may subsequently be composited together after each frequency in the minimal subset has 
been image processed. 



[0027] One of the programming problems in a parallel environment is that N as defined 
5 above may not equal the actual number of processors required to perform the task. Since 
the number of frequencies is determined by the ratio 

* = f (8) 

one need only adjust R so that N fits the initial number of processors exactly. After 

N 

10 N f is the smallest integer greater than — , compute 

N 



I log /, -log /q 

* = Vio N >* N + -l . 



(9) 



[0028] The new number, N, of frequencies is then assured to be N = N p * N fpe . This 

done, the number of frequencies per processor will be at least one and the number to 
15 process will divide evenly by the number of processors. As the process continues, the 

actual number of frequencies will decrease, and as a result, so will the required number of 
nodes. 



20 



[0029] To see how this process reduces the number of frequencies in a given migration, 
assume that the maximum offset is 6 km and the maximum migration depth is 10 km. If 
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the target frequency range is f 0 =5 and /j = 45 then, assuming a maximum angle of 
16.7 degrees 



+ 1 * 1.044 



the migration-frequency sequence looks something like 



5,5.22,5.45,5.69,5.94,6.2,6.47, . . . 



and 



TV = « 51. 

log 1.044 



(11) 



[0030] In contrast, with 8 seconds of data, the usual df would be — = 0.0625 and 

16 



40 

N freo =—^— + 1 = 641 (12) 
fm 0.0625 



The performance improvement is thus w 12.5 . 
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[0031] An even more dramatic result is possible by reducing the maximum migration 
depth from 10 km to 6 km. In this case, 



a x =VL25 « 1.12 (13) 



the migration-frequency set is 



5,5.6,6.27,7.025,7.87,... 



and 



A^J^i9. (14) 
logl.12 



641 

Now the performance gain is an astounding = 33.7 ! 



[0032] Figures 2, 3, 4, and 5 provide a series of images that demonstrate the limited 

frequency approach. They show that in some cases, one can achieve useable images with 

an extremely sparse data set. The first, Figure 2, is a fullband migration panel 201 of the 

SEG AA' data set using the full frequency band of 384 frequencies. The second panel 

labeled 301 in Figure 3 was produce with R= 58 and only 21 frequencies. The third panel 

401, in Figure 4, is the result with R=.292 for 73 frequencies. The fourth panel 501, 

Figure 5, is the result with R=.2 for 153 frequencies. Of these, only Figure 3 displays a 

questionable result. Depending on the processing needs, it might be acceptable as a 

starting point, but certainly not a final image. Figure 4 represents a good compromise 
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between speed and accuracy. Note that it was produce in 19% of the time it took for the 
fullband image in Figure 2. The final figure, Figure 5, is virtually identical to the 
broadband result in Figure 2. While it was produce in slightly less than half the time used 
for the full-band image it is not clear that the extra effort is worth the computational cost. 

[0033] The preceding figures demonstrate the validity of the frequency selection process 
and also show that, depending on the processing need, acceptable results can be achieved 
very efficiently. Limited- frequency methods (with minimal frequency subsets) provide 
for accurate imaging of large seismic data sets. If there is need to increase the density of 
frequencies defined by Equation (7), one simply starts with a new f Q , or ratio R, then 

remigrate with a new independent frequency set and composites or adds the results 
together. Note that the frequency density is completely controlled through the ratio of the 
offset wavenumber k h to the vertical wave number k 2 . As the incident angle approaches 

zero, fewer and fewer frequencies are required to image the same number of vertical and 
horizontal wavenumber. The fact that the number of frequencies in this implementation is 
on the order of at most 100, means that data distribution and resource allocation issues 
are minimized. The speed of the limited frequency approach means that routine 
application of full wave-equation methods is now a realistic, economical approach to 
high quality seismic imaging. 

[0034] Figure 6 illustrates the result of coupling intermediate frequency spline 
interpolation on a limited frequency downward continuation with R=.292. In this case, 
shot-domain migration 601 of the two-dimensional SEG AA' data set, the full frequency 

COR-1058 16 



band was computed from the limited frequency set by spline interpolation. The 
computational load in this case is roughly the same as that for Figure 4, but the results 
are at least as good. 

[0035] Figure 7 shows that the limited frequency approach works quite well on a depth 
slice 701 of real 3D data. This depth slice from a common-azimuth-prestack migration of 
a large Gulf of Mexico data set was obtained using an average of 100 frequency slices. 
The resulting runtime was approximately 1/9 that of a more traditional approach. 

[0036] The flow chart of Figure 8 illustrates a preferred embodiment of the method and 
system provided by the present invention. Seismic data are acquired 901 and transformed 
a frequency domain 903. The data may be gathered in any manner, for example common 
shot, common midpoint, common depth point, common image point, common offset, etc. 
A plurality of minimal frequency subsets are formed 905 comprising selected 
frequencies. The frequencies may be selected according to Equation 6 above so that 
more frequencies are selected for shallower image positions. The selected frequencies 
are imaged 907 to form frequency images. The selected frequencies may be individually 
weighted before or after imaging to assure a full and balanced bandwidth. The frequency 
images are composited 909 to form an intermediate or final composited migration image. 

[0037] The method and system of the present invention disclosed herein may be 
conveniently carried out by writing a computer program to carry out the steps described 
herein on a work s'ation or other conventional digital computer system of a type normally 
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used in the industry. The generation of such a program may be performed by those of 
ordinary skill in the art based on the processes described herein. Figure 9 illustrates a 
computer system comprising a central processing unit 1011, a display 1001, an input 
device 1021, and a plotter 1031. The computer program for carrying out the invention 
will normally reside on a storage media (not shown) associated with the central 
processing unit. The computer program may be transported on a CD-ROM or other 
storage media shown symbolically as storage medium 1041. 

[0038] The method and system of the present invention provides results that may be 
displayed or plotted with commercially available visualization software and computer 
peripherals. Such software and computer peripherals are well known to those of ordinary 
skill in the art. It should be appreciated that the results of the methods of the invention 
can be displayed, plotted and/or stored in various formats. 

[0039] While the invention has been described and illustrated herein by reference to 
embodiments in relation to the drawings attached hereto, various changes and further 
modifications, apart from those shown or suggested herein, may be made herein by those 
skilled in the art, without departing from the spirit of the invention, the scope of which is 
expressed in the following claims. 
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